# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
### Reducing Prejudice and Support for Religious 
### Nationalism Through Conversations on WhatsApp 

### Author: Rajeshwari Majumdar
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# File:     supp_robustness.R
# Purpose:  perform robustness checks described in Appendix D; 
#           produce summary stats, figures, and tables presented in 
#           all Appendix D subsections (unless stated otherwise below)
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

# Load data and custom functions 
rm(list=ls())
load("data/Completes_Merged_Anonymized.RData")
source("code/02_setup_functions.R")

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# D.1 Data quality checks ---- 
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

## Check how many flagged respondents
tabyl(d$maybe_duplicate_friend)

## Estimate main results (prejudice and religious nationalism) with flagged respondents removed
fn_reg_more(depvars = c("post.feelth_muslim","follow.feelth_muslim",
                          "twt_HinduNation_q1","twt_Hijab_q1","twt_HinduWeap_q1","twt_Terrorism_q1"),
              data = subset(d, maybe_duplicate_friend == 0)
)


# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# D.2 Balance checks ----
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

## see supp_balance.R


# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# D.3 Differential attrition ----
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

## see supp_participation.R


# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# D.4 Manipulation check ----
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

## Part 1: Characterize the problem 

### How many correct vs incorrect vs NA
tabyl(d$fq_religion_correct)

### Set NAs as incorrects
d$fq_religion_correct = ifelse(is.na(d$fq_religion_correct)==TRUE, 0, d$fq_religion_correct)
tabyl(d$fq_religion_correct)

### Hindu responses
dh = subset(d, religion=="Hindu")
t.test(dh$fq_religion_correct[dh$condition_pair=="HH"], dh$fq_religion_correct[dh$condition_pair=="HM"])

### Muslim responses
d %>% filter(religion == "Muslim", condition_pair == "HM") %>% summarise(mean(fq_religion_correct, na.rm = TRUE))

### Hindu responses by topic
d %>% filter(religion == "Hindu", condition_pair == "HM") %>% group_by(condition_topic) %>% summarise(mean(fq_religion_correct, na.rm = TRUE))

### Check how responses compare to the identification of other demographics
d %>%
  summarize(
    fq_religion_correct = mean(fq_religion_correct == 1, na.rm = TRUE),
    fq_age_correct      = mean(fq_age_correct == 1, na.rm = TRUE),
    fq_language_correct = mean(fq_language_correct == 1, na.rm = TRUE),
    fq_state_correct    = mean(fq_state_correct == 1, na.rm = TRUE)
  ) 

## Part 2: IV
## Table D.7 Effect of treatment on short-term feelings about Muslims (IV estimates) ----

### Subset to Hindu respondents 
dh = subset(d, religion=="Hindu")

### Convert treatment assignment variable to numeric
dh$condition_pair_num = as.numeric(as.factor(dh$condition_pair))

### Set up the treated variable: whether partner is perceived to be Muslim
dh$fq_religion_muslim = ifelse(dh$fq_religion_hm == "Muslim",1,0)
dh$fq_religion_muslim = ifelse(is.na(dh$fq_religion_muslim)==TRUE, 0, dh$fq_religion_muslim)
table(dh$condition_pair, dh$fq_religion_muslim)

### IV models for prejudice
iv.model = ivreg(post.feelth_muslim ~ fq_religion_muslim | condition_pair_num + prejudice_binary, data = dh)
iv.model.non = ivreg(post.feelth_muslim ~ fq_religion_muslim | condition_pair_num + prejudice_binary, 
                     data = subset(dh, condition_topic == "non"))
iv.model.gen = ivreg(post.feelth_muslim ~ fq_religion_muslim | condition_pair_num + prejudice_binary, 
                     data = subset(dh, condition_topic == "gen"))
iv.model.igr = ivreg(post.feelth_muslim ~ fq_religion_muslim | condition_pair_num + prejudice_binary, 
                     data = subset(dh, condition_topic == "igr"))

summary(iv.model, diagnostics = TRUE)

### Print
modelsummary(list(iv.model, iv.model.non, iv.model.gen, iv.model.igr),
             coef_map = c(
               "fq_religion_muslim" = "Identifies partner as Muslim",
               "(Intercept)" = "Constant"
             ),
             stars = c(`***` = 0.01, `**` = 0.05, `*` = 0.10),
             gof_map = "nobs"
)


# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# D.5 Multiple comparisons ----
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

## see results_religiousnatlm.R


# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# D.6 Alternate specifications ---- 
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

## see results_prejudice.R and results_religiousnatlm.R


# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# D.7 Weighted regressions ----
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

## Table D.13: Demographically weighted estimates ----
depvars = c("post.feelth_muslim","follow.feelth_muslim",
            "twt_HinduNation_q1","twt_Hijab_q1","twt_HinduWeap_q1","twt_Terrorism_q1")
depvars = which(names(d) %in% depvars)
mod = list()
for (i in 1:length(depvars)){
  d$depvar = d[,depvars[i]]
  fit = lm(depvar ~ condition_pair + prejudice_binary, 
           data = d %>% filter(religion=="Hindu"), 
           weights = survey_weight)
  mod[[i]] = fit
}
modelsummary(mod,
             coef_map = c(
               "condition_pairHM" = "Hindu-Muslim pair",
               "prejudice_binary" = "Baseline prejudice",
               "(Intercept)" = "Constant"
             ),
             stars = c(`***` = 0.01, `**` = 0.05, `*` = 0.10),
             gof_map = "nobs"
             )


# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# D.8 Placebo tests ----
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

## Feelings scales (Figure D.2) ----
# Estimate models corresponding to each group that feelings are elicited about (both one day later and two weeks later)
results = fn_extract_coefs(depvars = c("post.feelth_hindu","post.feelth_muslim","post.feelth_christian","post.feelth_brahmin","post.feelth_dalit","post.feelth_bjpvoter","post.feelth_incvoter","post.feelth_lgbtq","post.feelth_studentactiv","post.feelth_police",
                                         "follow.feelth_hindu","follow.feelth_muslim","follow.feelth_christian","follow.feelth_brahmin","follow.feelth_dalit","follow.feelth_bjpvoter","follow.feelth_incvoter","follow.feelth_lgbtq","follow.feelth_studentactiv","follow.feelth_police"),
                             depvar_labels = rep(c("Hindus","Muslims","Christians","Brahmins","Dalits","BJP voters","INC voters","LGBTQ people","Student activists","Police officers"),2),
                             timepoints = c(rep("One day post-treatment",10), rep("Two weeks post-treatment",10)),
                             controls = c("feelth_hindu","feelth_muslim","feelth_christian","feelth_brahmin","feelth_dalit","feelth_bjpvoter","feelth_incvoter","feelth_lgbtq","feelth_studentactiv","feelth_police"),
                             data = d)

# Check FDR correction 
results %>% mutate(p_fdr = p.adjust(p.value, method = "fdr"))

# Make coefficient plot
results$depvar_label = fct_rev(factor(results$depvar_label, levels = unique(results$depvar_label)))
ggplot(results %>% filter(timepoint == "One day post-treatment"), aes(x = depvar_label, y = estimate)) +
  geom_point() +
  geom_errorbar(aes(ymin = estimate - 1.96*std.error, ymax = estimate + 1.96*std.error), width = 0) +
  labs(x = "", y = expression(Coefficient~estimate~of~italic("Hindu-Muslim pair"))) + 
  scale_y_continuous(limits = c(-1,1), breaks = seq(-1,1,0.25)) +
  geom_hline(yintercept = 0, color = "gray", lty = 2) +
  coord_flip() + theme_classic(base_size = 16) 


## Politician statements (Figure D.3) ----
# Estimate models corresponding to each politician statement (both immediate and three weeks later)
results = fn_extract_coefs(depvars = c("twt_HinduNation_q1","twt_Hijab_q1","twt_NatlmText_q1","twt_LangHindi_q1","twt_ChinaBoycott_q1","twt_GBV_q1","twt_CovidUnity_q1","twt_CovidDoc_q1",
                                         "twt_HinduNationFollow_q1","twt_HijabFollow_q1","twt_HinduWeap_q1","twt_Terrorism_q1","twt_Housing_q1","twt_Education_q1","twt_LangAll_q1","twt_Astrology_q1","twt_RusAgainst_q1","twt_RusSupport_q1","twt_ChinaBoycottFollow_q1"),
                             depvar_labels = c("India is a Hindu nation","Hijab bans in schools","Nationalizing curricula","Hindi national language","Boycotting Chinese products","Gender-based violence","Unity during COVID","COVID healthcare workers",   
                                               "India is a Hindu nation","Hijab bans in schools","Violence for Hindu nation","Islam and terrorism","Religious segregation","Education infrastructure","Equality of languages","Astrology above science","India should condemn Russia","India should support Russia","Boycotting Chinese products"),
                             timepoints = c(rep("Statements shown on the last day of treatment",8), rep("Statements shown three weeks post-treatment",11)),
                             data = d)

# Check FDR correction 
results %>% group_by(timepoint) %>% mutate(p_fdr = p.adjust(p.value, method = "fdr")) %>% ungroup()

# Make coefficient plot
label_lookup = results %>% distinct(depvar, depvar_label)
ggplot(results, aes(x = reorder(depvar, -estimate), y = estimate)) +
  facet_wrap(~ timepoint, scales = "free_y") +
  geom_point() +
  geom_errorbar(aes(ymin = estimate - 1.96*std.error, ymax = estimate + 1.96*std.error), width = 0) +
  labs(x = "", y = expression(Coefficient~estimate~of~italic("Hindu-Muslim pair"))) + 
  scale_x_discrete(labels = function(x) {label_lookup$depvar_label[match(x, label_lookup$depvar)]}) +
  scale_y_continuous(limits = c(-0.3,0.3), breaks = seq(-0.3,0.3,0.1), labels = function(x) sprintf("%.1f", x)) +
  geom_hline(yintercept = 0, color = "gray", lty = 2) +
  coord_flip() + theme_classic(base_size = 16) + theme(strip.background = element_blank()) 


# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# D.9 Expected null effects ----
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

## Table D.14: Effect of treatment on opinions related to endogamy and affirmative action ----
fn_reg_more(depvars = c("opn_reservation_religion","comf_relationship_muslim","opn_interfaith_bad","opn_interfaith_ban"), add.controls = "prejudice_binary", data = d)
